program main

 implicit none

 call star ()

 stop
end

!************************************************************

subroutine star ()

 implicit none

 integer (kind=4),parameter :: n=2
 integer (kind=4), parameter :: l=64
 integer (kind=4) k
 external star_f
 real (kind=8) x0,x1,xmax,dx
 real (kind=8) u0(n),u1(n)
 real (kind=8) the0, the1
 real (kind=8) Rg,cp,cap

 x0=-0.2
 xmax=1.2
 dx=(xmax-x0)/(l-1)
 Rg=2./5.
 cp=1.
 cap=Rg/cp

 u0(1)=2.1
 u0(2)=150.
 the0=u0(1)*(2.1*150./(u0(1)*u0(2)))**cap

 do k=1,l
  write( *, '(2x,g14.6,2x,g14.6,2x,g14.6,2x,g16.8)') x0,u0(1),u0(2),the0 
  x1=x0+dx
  call rk4vec(x0,n,u0,dx,star_f,u1)
  the1=u1(1)*((2.1*150)/(u1(1)*u1(2)))**cap
 
  x0=x1
  u0(1:n)=u1(1:n)
  the0=the1

 end do
 
 return
end

!*************************************************************

subroutine star_f (x,n,u,uprime)

 implicit none
 
 integer (kind=4) n
 real (kind=8) x
 real (kind=8) u(n)
 real (kind=8) uprime(n)
 real (kind=8) g0
 real (kind=8) Rg
 real (kind=8) mr
 real (kind=8) m1
 real (kind=8) m2
 real (kind=8) m3
 real (kind=8) d
 real (kind=8) xtac
 real (kind=8) xtop
 
 Rg=2./5.
 g0=-1
 m1=2.
 m2=1.49
 m3=1.4
 d=0.05
 xtac=0.0
 xtop=0.9

 mr=m1 - (m1-m2)*0.5*(1.+erf((x-xtac)/(d))) - (m2-m3)*0.5*(1.+erf((x-xtop)/(d)))
 uprime(1)=g0/((mr+1)*Rg)
 uprime(2)=(u(2)/u(1))*(g0/Rg-uprime(1))
 
 return 
end

!**************************************************************

subroutine rk4vec (t0,m,u0,dt,f,u)

 implicit none

 integer (kind=4) m
 real(kind=8) dt
 external f
 real (kind=8) f0(m),f1(m),f2(m),f3(m)
 real (kind=8) t0,t1,t2,t3
 real (kind=8) u(m),u0(m),u1(m),u2(m),u3(m)

 call f (t0,m,u0,f0)
 
 t1=t0+dt/2.0D+00
 u1(1:m)=u0(1:m)+dt*f0(1:m)/2.0D+00
 call f (t1,m,u1,f1)

 t2=t0+dt/2.0D+00
 u2(1:m)=u0(1:m)+dt*f1(1:m)/2.0D+00
 call f (t2,m,u2,f2)

 t3=t0+dt
 u3(1:m)=u0(1:m)+dt*f2(1:m)
 call f (t3,m,u3,f3)

 u(1:m)=u0(1:m)+dt*(f0(1:m)+2.0D+00*f1(1:m)+2.0D+00*f2(1:m)+f3(1:m))/6.0D+00
 
 return
end

